{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## _*Particle hole transformation of FermionicOperator*_\n",
    "\n",
    "This notebook demonstrates carrying out a ParticleHole transformation on the FermionicOperator in Qiskit Chemistry. Here we use the FermionicOperator directly to demonstrate.\n",
    "\n",
    "Note: The Hamiltonian class that wraps this provides a means to use either full, or particle hole transformation. Under the covers it does what is shown here though.\n",
    "\n",
    "This notebook has been written to use the PYSCF chemistry driver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from qiskit import BasicAer\n",
    "\n",
    "from qiskit.aqua import QuantumInstance\n",
    "from qiskit.aqua.algorithms import VQE, NumPyMinimumEigensolver\n",
    "from qiskit.aqua.components.optimizers import L_BFGS_B\n",
    "from qiskit.circuit.library import TwoLocal\n",
    "\n",
    "from qiskit.chemistry import FermionicOperator\n",
    "from qiskit.chemistry.drivers import PySCFDriver, UnitsType"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll do this with H2 molecule and use the PySCF driver to create the integrals we need for the FermionicOperator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "driver = PySCFDriver(atom='H .0 .0 .0; H .0 .0 0.735', unit=UnitsType.ANGSTROM,\n",
    "                     charge=0, spin=0, basis='sto3g')\n",
    "molecule = driver.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create the FermionicOperator and use NumPyMinimumEigensolver with qubit operator we get from it via a jordan wigner mapping to compute the ground state energy. Here this is the electronic component of the total ground state energy (the total ground state energy would include the nuclear repulsion energy we can get from the molecule that comes from the driver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The exact ground state energy is: -1.8572750302023824\n",
      "The Hartree Fock Electron Energy is: -1.8369679912029842\n"
     ]
    }
   ],
   "source": [
    "ferOp = FermionicOperator(h1=molecule.one_body_integrals, h2=molecule.two_body_integrals)\n",
    "qubitOp_jw = ferOp.mapping(map_type='JORDAN_WIGNER', threshold=0.00000001)\n",
    "qubitOp_jw.chop(10**-10)\n",
    "\n",
    "# Using NumPyMinimumEigensolver to get the smallest eigenvalue\n",
    "exact_eigensolver = NumPyMinimumEigensolver(qubitOp_jw)\n",
    "ret = exact_eigensolver.run()\n",
    "\n",
    "print('The exact ground state energy is: {}'.format(ret.eigenvalue.real))\n",
    "print('The Hartree Fock Electron Energy is: {}'.format(molecule.hf_energy - molecule.nuclear_repulsion_energy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now the same as above but with ParticleHole transformation. This removes out energy from the FermionicOperator that is equivalent to the electronic part of the Hartree Fock Energy that we also computed above. The Hartree Fock energy also comes from the driver. To get the total electronic ground state energy we need to add the part we now compute with the part that was removed by the transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Energy shift is: 1.8369679912029837\n",
      "The exact ground state energy in PH basis is -0.020307038999395295\n",
      "The exact ground state energy in PH basis is -1.857275030202379 (with energy_shift)\n"
     ]
    }
   ],
   "source": [
    "# particle hole transformation\n",
    "newferOp, energy_shift = ferOp.particle_hole_transformation([molecule.num_alpha, molecule.num_beta])\n",
    "print('Energy shift is: {}'.format(energy_shift))\n",
    "newqubitOp_jw = newferOp.mapping(map_type='JORDAN_WIGNER', threshold=0.00000001)\n",
    "newqubitOp_jw.chop(10**-10)\n",
    "\n",
    "exact_eigensolver = NumPyMinimumEigensolver(newqubitOp_jw)\n",
    "ret = exact_eigensolver.run()\n",
    "\n",
    "print('The exact ground state energy in PH basis is {}'.format(ret.eigenvalue.real))\n",
    "print('The exact ground state energy in PH basis is {} (with energy_shift)'.format(ret.eigenvalue.real - energy_shift))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We run here using the quantum VQE algorithm to show the same result. The parameters printed are the optimal parameters of the variational form at the minimum energy, the ground state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimum value: -0.02030703892349\n",
      "Minimum value: -1.8572750301264738\n",
      "Parameters: [ 5.28462733  3.50885712  1.19186098 -1.31257921  5.55321399 -0.49555501\n",
      " -5.10169198 -5.4784798  -0.88598656  4.05399968  3.26356323  6.05812004\n",
      "  3.41375049  3.10439713 -1.10313698 -0.57661347  2.65713956 -1.90680973\n",
      "  6.10125861  0.49985016  4.83574309 -2.15677133 -2.96541528 -2.28211722]\n"
     ]
    }
   ],
   "source": [
    "# setup VQE \n",
    "# setup optimizer, use L_BFGS_B optimizer for example\n",
    "lbfgs = L_BFGS_B(maxfun=1000, factr=10, iprint=10)\n",
    "\n",
    "# setup variational form generator (generate trial circuits for VQE)\n",
    "var_form = TwoLocal(newqubitOp_jw.num_qubits, 'ry', 'cz', reps=5, entanglement=[[0, 1], [1, 2], [2, 3]])\n",
    "\n",
    "# setup VQE with operator, variational form, and optimizer\n",
    "vqe_algorithm = VQE(newqubitOp_jw, var_form, lbfgs)\n",
    "\n",
    "backend = BasicAer.get_backend('statevector_simulator')\n",
    "quantum_instance = QuantumInstance(backend)\n",
    "\n",
    "results = vqe_algorithm.run(quantum_instance)\n",
    "print(\"Minimum value: {}\".format(results.eigenvalue.real))\n",
    "print(\"Minimum value: {}\".format(results.eigenvalue.real - energy_shift))\n",
    "print(\"Parameters: {}\".format(results.optimal_point))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
